Treatment order is not a significant predictor of wingbeat frequency in the model with percent load (though it is very close, p ~ 0.07)
percent load and treatment suggest very similar models in terms of wingbeat frequency
Treatment order is not a strong predictor of stroke amplitude in the model with percent load
wingbeat freq = size Met rate = mass carrying amplitude = percent loading
Susie’s (hypotheses): lm( frequency ~ ITspan + trt order)
lm( metabolism ~ total mass)
lm( stroke amplitude ~ %load + trt) (those are probably the same thing)
How heavy the bee is loaded overall, affects how much the amplitude and frequency change
Questions to look into: ** A really heavy loaded bee will have a smaller change in metabolism vs. light bee ** Suggests a non-linear relationship or random slopes that depend on bee size — the cost for additional weight goes down as you get heavier loaded. — not the case for the change in amplitude (change in load vs. change in amplitude is linear) — not the case for metabolism or wingbeat frequency.
** Hypotheses: 1. Wingbeat freq moves around for a number of reasons (tired, cold, added mass) 2. Stroke amplitude rescues wingbeat frequency 3. As you get really close to your max, adding more weight doesn’t seem to cost as much.
Analyses of respirometry data:
Install required packages and read in data Define custom function for evaluating VIF with multilevel models
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS")
ipak(packages)
## ggplot2 car lme4 gsheet MASS
## TRUE TRUE TRUE TRUE TRUE
theme_set(theme_bw())
# read in data -- google sheet called "Bumble mumble grumble"
bdta <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1GUUoFq41Ep3sJNFXiyMcp7fZhYmQCzTcLnhVXdr4WRo/edit?usp=sharing")
(used later)
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}
# calculate percent load
bdta$percLoad <- with(bdta, (Mstarved + load) / Mstarved)
# calculate total mass
bdta$totalMass <- with(bdta, Mstarved + load)
summary(bdta)
## Bee.ID Treatment.order Treatment Mstarved
## Length:60 Min. :1.0 Length:60 Min. :0.0767
## Class :character 1st Qu.:1.0 Class :character 1st Qu.:0.1213
## Mode :character Median :1.5 Mode :character Median :0.1378
## Mean :1.5 Mean :0.1411
## 3rd Qu.:2.0 3rd Qu.:0.1670
## Max. :2.0 Max. :0.1909
## load IT.Span single.wing.area av.resp..CO2.mL.hr.
## Min. :0.01650 Min. :3.860 Min. :20.03 Min. : 5.507
## 1st Qu.:0.04713 1st Qu.:4.400 1st Qu.:27.93 1st Qu.: 8.935
## Median :0.07840 Median :4.620 Median :30.67 Median :10.666
## Mean :0.08332 Mean :4.594 Mean :30.70 Mean :10.600
## 3rd Qu.:0.11800 3rd Qu.:4.870 3rd Qu.:33.75 3rd Qu.:12.297
## Max. :0.17460 Max. :5.180 Max. :39.05 Max. :16.475
## wbf.aud. stroke.amplitude wing.velocity percLoad
## Min. :163.9 Min. : 96.88 Min. :1.208 Min. :1.143
## 1st Qu.:173.0 1st Qu.:114.16 1st Qu.:1.351 1st Qu.:1.321
## Median :179.1 Median :125.96 Median :1.427 Median :1.576
## Mean :178.9 Mean :123.31 Mean :1.465 Mean :1.600
## 3rd Qu.:185.3 3rd Qu.:133.61 3rd Qu.:1.588 3rd Qu.:1.863
## Max. :194.7 Max. :144.78 Max. :1.843 Max. :2.180
## totalMass
## Min. :0.1062
## 1st Qu.:0.1808
## Median :0.2191
## Mean :0.2245
## 3rd Qu.:0.2646
## Max. :0.3550
ggplot(bdta, aes(x = percLoad, fill = Treatment)) +
geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# seems like percent load and trt are very related
# visualize % load vs total mass
ggplot(bdta, aes(x = percLoad, y = Mstarved + load)) +
geom_point()
# seem to be strongly associated
# see number of observations per bee
table(bdta$Bee.ID) # each bee has two observations
##
## E01 E03 E05 E09 E10 E11 E12 E16 E20 E24 E26 E28 E30 E31 E32 E33 E35 E36
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## E37 E39 E41 E42 E44 E45 E49 E50 E53 E54 E56 E59
## 2 2 2 2 2 2 2 2 2 2 2 2
# get number of unque bees
length(table(bdta$Bee.ID))
## [1] 30
# get treatment orders
trtOrders_loadedSecond <- sapply(X = unique(bdta$Bee.ID), FUN = function(x){
tmp = bdta[bdta$Bee.ID == x, c("Treatment.order", "Treatment")]
if("2_L" %in% paste(tmp[,1], tmp[,2], sep = "_")){
loadedSecond = TRUE
}
else loadedSecond = FALSE
return(loadedSecond)
})
table(trtOrders_loadedSecond)
## trtOrders_loadedSecond
## FALSE
## 30
# visualize scatterplot
car::scatterplotMatrix(bdta[, 4:13])
# VIF is high among the bee size predictors
# note: this function is almost the same as car::vif
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + Mstarved + IT.Span + totalMass + single.wing.area + percLoad + (1|Bee.ID), data = bdta))
## TreatmentUL Treatment.order Mstarved IT.Span
## 16.591645 1.270430 23.112962 6.658432
## totalMass single.wing.area percLoad
## 32.931859 10.903672 38.415618
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area")])
# principle components
aa = prcomp(bdta[, c("Mstarved", "IT.Span", "single.wing.area")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
## PC1 PC2 PC3
## Standard deviation 1.6852 0.33843 0.21323
## Proportion of Variance 0.9467 0.03818 0.01516
## Cumulative Proportion 0.9467 0.98484 1.00000
biplot(aa) # shows that all three size measurement are correlated
# note, I changed the signs of the predictions so that higher PC1 values
# correspond to bigger bees
p1 = -predict(aa)[,1]
# add PC1 scores to dataset
bdta$size_pca1 = p1
# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area", "size_pca1", "percLoad")])
# check VIF one more time
# VIF is high among the bee size predictors
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + size_pca1 + percLoad + totalMass + (1|Bee.ID), data = bdta))
## TreatmentUL Treatment.order size_pca1 percLoad
## 16.643643 1.258990 7.448179 30.631628
## totalMass
## 26.790114
# remove treatment
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pca1 + percLoad + totalMass + (1|Bee.ID), data = bdta))
## Treatment.order size_pca1 percLoad totalMass
## 1.014834 7.698050 20.910198 25.813622
# looks like we cannot have total mass and percentLoad in the same model, according to VIF
# remove total mass
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pca1 + percLoad + (1|Bee.ID), data = bdta))
## Treatment.order size_pca1 percLoad
## 1.003811 1.004814 1.008624
# a few house-keeping issues to make models easier to read
# convert trt order to factor
bdta$Treatment.order <- as.factor(as.character(bdta$Treatment.order))
# change reference level to unloaded
bdta$Treatment <- factor(bdta$Treatment, levels = c("UL", "L"))
# make a full model with all two-way interactions
m1 <- lmer(av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 196.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5230 -0.4541 -0.0409 0.4813 2.2994
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.2498 1.1179
## Residual 0.7711 0.8781
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.0096 1.3291 2.264
## size_pca1 0.6926 0.4455 1.555
## Treatment.order2 1.9553 2.0573 0.950
## percLoad 4.8854 0.8030 6.084
## size_pca1:Treatment.order2 -0.0830 0.1451 -0.572
## size_pca1:percLoad 0.3533 0.2491 1.418
## Treatment.order2:percLoad -1.4903 1.2730 -1.171
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.223
## Trtmnt.rdr2 -0.860 0.246
## percLoad -0.980 0.204 0.866
## sz_pc1:Tr.2 0.297 -0.379 -0.273 -0.297
## sz_pc1:prcL 0.148 -0.933 -0.184 -0.125 0.228
## Trtmnt.r2:L 0.853 -0.233 -0.994 -0.868 0.271 0.169
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad)
anova(m1, m2.0) # likelihood ratio test for interaction of treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.0: (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 208.02 224.78 -96.011 192.02
## m1 9 208.55 227.40 -95.275 190.55 1.472 1 0.225
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 199.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61501 -0.46662 -0.04214 0.53248 2.31063
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.3003 1.1403
## Residual 0.7595 0.8715
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.33950 0.69115 6.279
## size_pca1 0.57136 0.43113 1.325
## Treatment.order2 -0.43829 0.22783 -1.924
## percLoad 4.06743 0.39578 10.277
## size_pca1:Treatment.order2 -0.03698 0.13862 -0.267
## size_pca1:percLoad 0.40235 0.24379 1.650
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.048
## Trtmnt.rdr2 -0.219 0.137
## percLoad -0.924 0.003 0.055
## sz_pc1:Tr.2 0.131 -0.337 -0.036 -0.131
## sz_pc1:prcL 0.007 -0.930 -0.144 0.045 0.192
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.1: (1 | Bee.ID) + size_pca1:percLoad
## m2.0: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.0: (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 206.10 220.76 -96.051 192.10
## m2.0 8 208.02 224.78 -96.011 192.02 0.0814 1 0.7754
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## (1 | Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 197.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.68107 -0.47084 -0.05908 0.51796 2.31287
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.3135 1.1461
## Residual 0.7334 0.8564
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.3659 0.6752 6.466
## size_pca1 0.5329 0.4000 1.332
## Treatment.order2 -0.4405 0.2237 -1.969
## percLoad 4.0521 0.3858 10.502
## size_pca1:percLoad 0.4146 0.2352 1.763
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.004
## Trtmnt.rdr2 -0.216 0.132
## percLoad -0.921 -0.044 0.051
## sz_pc1:prcL -0.019 -0.935 -0.140 0.072
m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.2: (1 | Bee.ID)
## m2.1: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.1: (1 | Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 207.38 219.94 -97.688 195.38
## m2.1 7 206.10 220.76 -96.051 192.10 3.2723 1 0.07046 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 199.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.44974 -0.62601 -0.00733 0.59569 2.08255
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.296 1.138
## Residual 0.785 0.886
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.3843 0.6950 6.309
## size_pca1 1.1922 0.1423 8.378
## Treatment.order2 -0.3854 0.2292 -1.681
## percLoad 4.0055 0.3976 10.074
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.064
## Trtmnt.rdr2 -0.221 0.004
## percLoad -0.925 0.069 0.062
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order (different than model with treatment)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 5 208.25 218.73 -99.127 198.25
## m2 6 207.38 219.94 -97.688 195.38 2.8789 1 0.08975 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for size
m4 <- update(m3, .~. - size_pca1)
anova(m3, m4) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: av.resp..CO2.mL.hr. ~ percLoad + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 4 244.02 252.40 -118.011 236.02
## m3 5 208.25 218.73 -99.127 198.25 37.767 1 7.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for Treatment (load)
m5 <- update(m3, .~. - percLoad)
anova(m3, m5) # keep percent load
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: av.resp..CO2.mL.hr. ~ size_pca1 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 255.27 263.64 -123.634 247.27
## m3 5 208.25 218.73 -99.127 198.25 49.013 1 2.543e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 201.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.33113 -0.59393 0.02882 0.58497 2.24345
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.2695 1.1267
## Residual 0.8349 0.9137
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.1216 0.6955 5.926
## size_pca1 1.1932 0.1423 8.388
## percLoad 4.0493 0.4087 9.907
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pca1 -0.067
## percLoad -0.940 0.071
# write output
summary(m3)$coefficients
## Estimate Std. Error t value
## (Intercept) 4.121569 0.6955406 5.925706
## size_pca1 1.193243 0.1422611 8.387692
## percLoad 4.049269 0.4087169 9.907270
write.csv( summary(m3)$coefficients, file = "RespCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good
# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good
# make a full model with all two-way interactions
m1 <- lmer(wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1 |
## Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 349.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.80184 -0.60216 -0.06242 0.61714 1.64271
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 10.18 3.190
## Residual 20.62 4.541
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 173.7585 5.7634 30.149
## size_pca1 -4.8135 2.1942 -2.194
## Treatment.order2 -5.8644 8.3595 -0.702
## percLoad 4.3885 3.4953 1.256
## size_pca1:Treatment.order2 0.9035 0.7397 1.221
## size_pca1:percLoad 0.9858 1.2643 0.780
## Treatment.order2:percLoad 1.4493 5.1496 0.281
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.198
## Trtmnt.rdr2 -0.818 0.226
## percLoad -0.984 0.172 0.814
## sz_pc1:Tr.2 0.254 -0.376 -0.223 -0.252
## sz_pc1:prcL 0.127 -0.958 -0.171 -0.098 0.216
## Trtmnt.r2:L 0.807 -0.209 -0.990 -0.818 0.221 0.152
### LRT's for interactions
m2.0 <- update(m1, .~. - size_pca1:percLoad )
anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.0: size_pca1:Treatment.order + Treatment.order:percLoad
## m1: wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1 |
## m1: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 381.69 398.44 -182.84 365.69
## m1 9 382.99 401.83 -182.49 364.99 0.7005 1 0.4026
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## size_pca1:Treatment.order + Treatment.order:percLoad
## Data: bdta
##
## REML criterion at convergence: 352.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72938 -0.54168 0.01881 0.54877 1.67034
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 10.08 3.175
## Residual 20.49 4.526
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 173.1944 5.6964 30.404
## size_pca1 -3.1748 0.6277 -5.058
## Treatment.order2 -4.7581 8.2056 -0.580
## percLoad 4.6533 3.4658 1.343
## size_pca1:Treatment.order2 0.7793 0.7199 1.082
## Treatment.order2:percLoad 0.8413 5.0701 0.166
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.268
## Trtmnt.rdr2 -0.815 0.221
## percLoad -0.984 0.272 0.813
## sz_pc1:Tr.2 0.234 -0.606 -0.194 -0.238
## Trtmnt.r2:L 0.803 -0.222 -0.990 -0.816 0.194
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.1: Treatment.order:percLoad
## m2.0: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.0: size_pca1:Treatment.order + Treatment.order:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 380.97 395.64 -183.49 366.97
## m2.0 8 381.69 398.44 -182.84 365.69 1.2893 1 0.2562
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## Treatment.order:percLoad
## Data: bdta
##
## REML criterion at convergence: 354.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.60845 -0.62485 0.06463 0.64370 1.69662
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 10.35 3.218
## Residual 20.39 4.515
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.6958 5.5454 30.962
## size_pca1 -2.7623 0.5021 -5.501
## Treatment.order2 -2.9560 8.0745 -0.366
## percLoad 5.5797 3.3706 1.655
## Treatment.order2:percLoad -0.2770 4.9889 -0.056
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.163
## Trtmnt.rdr2 -0.808 0.133
## percLoad -0.983 0.165 0.806
## Trtmnt.r2:L 0.796 -0.133 -0.989 -0.809
m2.2 <- update(m2.1, .~. - Treatment.order:percLoad)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## m2.1: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.1: Treatment.order:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 378.98 391.54 -183.49 366.98
## m2.1 7 380.97 395.64 -183.49 366.97 0.0018 1 0.9666
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 360
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64247 -0.64283 0.06575 0.64332 1.72090
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 9.874 3.142
## Residual 20.221 4.497
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.9557 3.3390 51.50
## size_pca1 -2.7663 0.4908 -5.64
## Treatment.order2 -3.3999 1.1632 -2.92
## percLoad 5.4191 1.9692 2.75
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.095
## Trtmnt.rdr2 -0.231 0.006
## percLoad -0.954 0.099 0.060
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi")
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wbf.aud. ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 5 384.96 395.43 -187.48 374.96
## m2 6 378.98 391.54 -183.49 366.98 7.9842 1 0.004719 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for size
m4 <- update(m2, .~. - size_pca1)
anova(m2, m4)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: wbf.aud. ~ Treatment.order + percLoad + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 399.76 410.23 -194.88 389.76
## m2 6 378.98 391.54 -183.49 366.98 22.783 1 1.813e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for Treatment (load)
m5 <- update(m2, .~. - percLoad)
anova(m2, m5)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: wbf.aud. ~ size_pca1 + Treatment.order + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 5 383.97 394.45 -186.99 373.97
## m2 6 378.98 391.54 -183.49 366.98 6.9967 1 0.008166 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 360
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64247 -0.64283 0.06575 0.64332 1.72090
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 9.874 3.142
## Residual 20.221 4.497
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.9557 3.3390 51.50
## size_pca1 -2.7663 0.4908 -5.64
## Treatment.order2 -3.3999 1.1632 -2.92
## percLoad 5.4191 1.9692 2.75
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.095
## Trtmnt.rdr2 -0.231 0.006
## percLoad -0.954 0.099 0.060
# write output
summary(m2)$coefficients
## Estimate Std. Error t value
## (Intercept) 171.955651 3.3390335 51.498630
## size_pca1 -2.766273 0.4908262 -5.635952
## Treatment.order2 -3.399880 1.1631609 -2.922966
## percLoad 5.419093 1.9691779 2.751957
write.csv( summary(m2)$coefficients, file = "FreqCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m2), main = "")
qqline(resid(m2)) # good
# residual plot
plot(fitted(m2), resid(m2), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m2)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m2)$Bee.ID[[1]]) # looks good
# make a full model with all two-way interactions
m1 <- lmer(stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 359.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89733 -0.52961 -0.01694 0.40969 1.99945
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.12 3.758
## Residual 23.72 4.870
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 61.9443 6.3305 9.785
## size_pca1 6.0503 2.3672 2.556
## Treatment.order2 7.4111 9.2790 0.799
## percLoad 38.4339 3.8380 10.014
## size_pca1:Treatment.order2 -1.1440 0.7947 -1.440
## size_pca1:percLoad -2.6118 1.3601 -1.920
## Treatment.order2:percLoad -5.0714 5.7204 -0.887
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.201
## Trtmnt.rdr2 -0.825 0.229
## percLoad -0.984 0.176 0.822
## sz_pc1:Tr.2 0.259 -0.377 -0.230 -0.258
## sz_pc1:prcL 0.129 -0.956 -0.172 -0.102 0.217
## Trtmnt.r2:L 0.814 -0.212 -0.991 -0.826 0.227 0.154
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad )
anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 393.14 409.89 -188.57 377.14
## m1 9 394.31 413.16 -188.16 376.31 0.8251 1 0.3637
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 365.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89468 -0.52361 -0.02085 0.43365 2.07383
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.75 3.841
## Residual 23.18 4.814
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.4391 3.6420 18.242
## size_pca1 5.5979 2.2913 2.443
## Treatment.order2 -0.7365 1.2581 -0.585
## percLoad 35.6715 2.1438 16.639
## size_pca1:Treatment.order2 -0.9857 0.7651 -1.288
## size_pca1:percLoad -2.4201 1.3296 -1.820
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.050
## Trtmnt.rdr2 -0.227 0.138
## percLoad -0.950 0.002 0.054
## sz_pc1:Tr.2 0.132 -0.345 -0.035 -0.128
## sz_pc1:prcL 0.006 -0.955 -0.142 0.046 0.189
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## m2.0: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 392.99 407.65 -189.49 378.99
## m2.0 8 393.14 409.89 -188.57 377.14 1.8488 1 0.1739
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 368.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92050 -0.50058 -0.06054 0.48999 2.08830
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.24 3.773
## Residual 23.88 4.887
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.1284 3.6558 18.362
## size_pca1 4.5889 2.1787 2.106
## Treatment.order2 -0.7942 1.2763 -0.622
## percLoad 35.2724 2.1552 16.367
## size_pca1:percLoad -2.1023 1.3241 -1.588
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.005
## Trtmnt.rdr2 -0.225 0.135
## percLoad -0.950 -0.046 0.050
## sz_pc1:prcL -0.019 -0.966 -0.138 0.072
m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.2: Bee.ID)
## m2.1: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 393.68 406.24 -190.84 381.68
## m2.1 7 392.99 407.65 -189.49 378.99 2.6909 1 0.1009
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 373.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87059 -0.50457 -0.03627 0.58077 2.02395
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.31 3.782
## Residual 24.73 4.973
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.0491 3.7152 18.047
## size_pca1 1.2463 0.5668 2.199
## Treatment.order2 -1.0742 1.2864 -0.835
## percLoad 35.5004 2.1858 16.241
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.091
## Trtmnt.rdr2 -0.230 0.006
## percLoad -0.952 0.096 0.060
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## size
m3 <- update(m2, .~. - size_pca1)
anova(m2, m3, test = "Chi") # keep size_pca1
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: stroke.amplitude ~ Treatment.order + percLoad + (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 5 396.50 406.97 -193.25 386.50
## m2 6 393.68 406.24 -190.84 381.68 4.8203 1 0.02813 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ Treatment.order + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 379
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79024 -0.53442 0.00173 0.52524 1.85491
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.63 4.199
## Residual 24.87 4.987
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.501 3.741 18.04
## Treatment.order2 -1.084 1.290 -0.84
## percLoad 35.221 2.193 16.06
##
## Correlation of Fixed Effects:
## (Intr) Trtm.2
## Trtmnt.rdr2 -0.229
## percLoad -0.948 0.060
# trt order
m4 <- update(m2, .~. - Treatment.order)
anova(m2, m4) # drop treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 392.42 402.89 -191.21 382.42
## m2 6 393.68 406.24 -190.84 381.68 0.7374 1 0.3905
# LRT for Treatment (load)
m5 <- update(m4, .~. - percLoad)
anova(m4, m5)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: stroke.amplitude ~ size_pca1 + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 478.07 486.45 -235.04 470.07
## m4 5 392.42 402.89 -191.21 382.42 87.657 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 376.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96941 -0.48117 -0.05909 0.48300 1.94459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.53 3.811
## Residual 24.42 4.941
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.3057 3.5966 18.436
## size_pca1 1.2495 0.5675 2.202
## percLoad 35.6293 2.1693 16.424
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pca1 -0.091
## percLoad -0.965 0.095
# write output
summary(m4)$coefficients
## Estimate Std. Error t value
## (Intercept) 66.305690 3.5965822 18.435750
## size_pca1 1.249463 0.5674547 2.201872
## percLoad 35.629329 2.1693080 16.424283
write.csv( summary(m4)$coefficients, file = "AmpCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # ok
# residual plot
plot(fitted(m4), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good
# make a full model with all two-way interactions
m1 <- lmer(wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -78.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00651 -0.41257 -0.06598 0.54758 2.20129
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003893 0.06239
## Residual 0.005935 0.07704
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.14983 0.10145 21.190
## size_pca1 -0.12816 0.03757 -3.411
## Treatment.order2 -0.11249 0.14953 -0.752
## percLoad -0.41933 0.06150 -6.819
## size_pca1:Treatment.order2 0.01603 0.01258 1.274
## size_pca1:percLoad 0.05043 0.02155 2.340
## Treatment.order2:percLoad 0.05830 0.09222 0.632
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.203
## Trtmnt.rdr2 -0.828 0.230
## percLoad -0.984 0.178 0.826
## sz_pc1:Tr.2 0.262 -0.377 -0.233 -0.261
## sz_pc1:prcL 0.131 -0.954 -0.173 -0.103 0.218
## Trtmnt.r2:L 0.818 -0.214 -0.991 -0.830 0.231 0.155
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad )
anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 -103.25 -86.492 59.624 -119.25
## m1 9 -101.66 -82.811 59.830 -119.66 0.4129 1 0.5205
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -80.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07164 -0.43946 -0.00926 0.53663 2.21665
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003959 0.06292
## Residual 0.005800 0.07616
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.09780 0.05777 36.31
## size_pca1 -0.12306 0.03633 -3.39
## Treatment.order2 -0.01883 0.01990 -0.95
## percLoad -0.38734 0.03397 -11.40
## size_pca1:Treatment.order2 0.01421 0.01210 1.17
## size_pca1:percLoad 0.04829 0.02106 2.29
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.050
## Trtmnt.rdr2 -0.227 0.138
## percLoad -0.949 0.002 0.054
## sz_pc1:Tr.2 0.132 -0.345 -0.035 -0.129
## sz_pc1:prcL 0.006 -0.954 -0.142 0.046 0.190
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## m2.0: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 -103.70 -89.044 58.852 -117.70
## m2.0 8 -103.25 -86.492 59.624 -119.25 1.5428 1 0.2142
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -86.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.08227 -0.51145 -0.04384 0.65437 2.23563
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003885 0.06233
## Residual 0.005902 0.07682
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.08849 0.05768 36.21
## size_pca1 -0.10840 0.03436 -3.16
## Treatment.order2 -0.01800 0.02006 -0.90
## percLoad -0.38198 0.03395 -11.25
## size_pca1:percLoad 0.04363 0.02084 2.09
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.005
## Trtmnt.rdr2 -0.224 0.135
## percLoad -0.948 -0.046 0.050
## sz_pc1:prcL -0.019 -0.965 -0.138 0.072
m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## don't drop size:percLoad interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.2: Bee.ID)
## m2.1: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 -101.16 -88.596 56.581 -113.16
## m2.1 7 -103.70 -89.044 58.852 -117.70 4.5421 1 0.03307 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# renaming model to simplify later typing
m2 <- m2.1
##### LRTs for main effects
# trt order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## m2: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2: Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 -104.83 -92.261 58.413 -116.83
## m2 7 -103.70 -89.044 58.852 -117.70 0.8775 1 0.3489
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -91.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18903 -0.53485 -0.00798 0.63629 2.12781
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003898 0.06243
## Residual 0.005865 0.07658
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.07698 0.05606 37.05
## size_pca1 -0.10424 0.03395 -3.07
## percLoad -0.38053 0.03381 -11.25
## size_pca1:percLoad 0.04104 0.02058 1.99
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 percLd
## size_pca1 0.026
## percLoad -0.963 -0.053
## sz_pc1:prcL -0.052 -0.964 0.080
# summarize final model for paper
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -91.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18903 -0.53485 -0.00798 0.63629 2.12781
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003898 0.06243
## Residual 0.005865 0.07658
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.07698 0.05606 37.05
## size_pca1 -0.10424 0.03395 -3.07
## percLoad -0.38053 0.03381 -11.25
## size_pca1:percLoad 0.04104 0.02058 1.99
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 percLd
## size_pca1 0.026
## percLoad -0.963 -0.053
## sz_pc1:prcL -0.052 -0.964 0.080
# write output
summary(m3)$coefficients
## Estimate Std. Error t value
## (Intercept) 2.07697845 0.05606361 37.046819
## size_pca1 -0.10424096 0.03394778 -3.070627
## percLoad -0.38052624 0.03381059 -11.254646
## size_pca1:percLoad 0.04104281 0.02058223 1.994090
write.csv( summary(m3)$coefficients, file = "VelocCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good
# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good